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Abstract. We propose a new method for determining the 
stochastic or ordered nature of trajectories in non-integrable 
Hamiltonian dynamical systems. The method consists of con- 
structing a time-series from the divergence of nearby trajecto- 
ries and then performing a power spectrum analysis of the se- 
ries. Ordered trajectories present a spectrum that consists of a 
few spikes while the spectrum of stochastic trajectories is con- 
tinuous. A test of the method with three different systems, a 
2-D mapping as well as a 2-D and a 3-D Hamiltonian, shows 
that the method is fast and efficient, even in the case of sticky 
trajectories. The method is also applied to the motion of aster- 
oids in the Solar System. 
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1. Introduction 

The problem of distinguishing a chaotic from an ordered trajec- 
tory in a non-integrable Hamiltonian system has been a topic 
of active investigation since the pioneering work of Henon and 
Heiles (1964). Initially, when the study was restricted to 2-D 
systems, the work was done through surface of section plots. 
Later, when systems with more than 2-D were considered, the 
method of choice was the calculation of the Lyapunov Char- 
acteristic Numbers (LCNs) (Benettin et al. 1976, Froeschle 
1984). Unfortunately both the above methods suffer from the 
same drawback, namely they are not able to distinguish easily 
an ordered from a "sticky" chaotic trajectory. Various methods 
have been devised since then to address the above problem, 
namely the distinction of an ordered from a chaotic trajectory 
using a relatively short-time trajectory segment. 

These methods generally fall into two main classes: those 
that use frequency or correlation analysis of a time-series, con- 
structed by the values of generalized co-ordinates (or functions 
of them), and those that use the geodesic divergence of initially 
nearby trajectories. In the first class belong the "old" method 
of the rotation number (Contopoulos 1966), the frequency map 
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analysis developed by Laskar (Laskar et al. 1992, Laskar 1993) 
and the power spectrum analysis of quasi-integrals developed 
by Voyatzis & Ichtiaroglou (1992). In the second belong the 
probability density function analysis of stretching numbers de- 
veloped by Froeschle et al. (1993) and Voglis & Contopoulos 
(1994) and the Fast Lyapunov Indicators method developed by 
Froeschle et al. (1997). Each one of the above methods has its 
own advantages and weaknesses; in particular some of them are 
more suitable to test large sets of trajectories rather than single 
ones, some are more efficient for 2-D systems rather than for 
N(>2)-D systems and some perform better for mappings rather 
than flows. 

In 1997, Contopoulos & Voglis introduced a new method 
for distinguishing chaotic from ordered trajectories, which does 
not belong to any of the above mentioned two classes but, 
instead, may be classified as "mixed". This new method is 
based on the analysis of the probability density of helicity 
and twist angles. Voglis & Efthymiopoulos (1998) and subse- 
quently Froeschle & Lega (1998) showed that the twist angles 
method is very efficient in testing phase space regions, at least 
in cases of 2-D systems, where the twist angles can be eas- 
ily calculated. More recently Voglis et al. (1998, 1999) pro- 
posed two new and very efficient methods, namely the method 
of "Dynamical Spectral Distance" (DSD), which is particularly 
suitable for the characterization of single trajectories in 4-D 
maps, and the method of "Rotational Tori Recognizer" (RO- 
TOR) which is very efficient for testing wide areas of the phase 
space of 2-D maps. 

In the present paper we are introducing a new "mixed" 
method, which we show that it is at least as sensitive as the 
other methods in the literature, may be applied in a straightfor- 
ward way to dynamical systems with more than two degrees of 
freedom and is equally efficient for single trajectories as well 
as large sets of them. The method consists in analyzing a time 
series constructed by the values of the geodesic deviation of 
nearby trajectories recorded at a properly selected frequency. It 
should be noted that a method based on a similar technique has 
been proposed by Lohinger and Froeschle (1993). 

The paper is organized as follows. Section 2 describes the 
basic features of the method, which, in Section 3, is tested upon 
three dynamical systems of different types. A comparison of 
the results of our method to those derived by various other 
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methods is made in Section 4. Section 5 treats the application 
of the proposed method to one of the most important problems 
of solar system dynamics, the motion of asteroids. Finally in 
Section 6 we present our conclusions. 

2. The method 

In order to decide on the nature of a trajectory (chaotic or not), 
we work as follows. We integrate numerically the "main" tra- 
jectory together with a nearby one, which at a time t = starts 
at an infinitesimal distance in phase space, do, from the main, 
and we calculate their distance, d\, at a time t — to + At. Let 
us denote by q the logarithm of the ratio of the two distances, 
q = \n(di/do). We then renormalize the nearby trajectory, so 
as to start from a new position in phase space, which is at dis- 
tance d from the main trajectory in the direction of d\. The 
trajectories are followed once again for a time interval At and 
a new q is calculated. After following the trajectories for a time 
interval T = N At, we have constructed a time series consist- 
ing of the consecutive q's 



q(t) = In 



do 



or q k 



In 



do 



(1) 



An ordered trajectory of a N-D conservative dynamical sys- 
tem will lie on an invariant torus, i.e. a N-D manifold of the 
2N-D phase space. Any such trajectory is, in general, quasi- 
periodic and covers densely the invariant set. If a nearby tra- 
jectory is started at an infinitesimal distance from the previous 
one, this too would, in general, lie on an invariant torus, so that 
the q(t) time series should behave in a quasi-periodic manner. 
On the other hand, a chaotic trajectory visits different regions 
of phase space in a stochastic manner and q(t) should also be 
"random". 

Following the above considerations we calculate the power 
spectrum P(f) of the q(t) time series. We first calculate the 
discrete Fourier transform of the q k multiplied by a window 
function w k 



N-l 



Qj = ]T q k w k e 2 ^ k / N j = 0..N - 1 
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Then the power spectrum, P(f), is defined in M — N/2 
frequencies as 
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and f c = f N / 2 is the Nyquist frequency defined as 

f c =—. 
Jc 2 At 
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Fig. 1. The successive points of the "sticky" trajectory (starting 
at Jo = 7T, 8 — 1.538 ir). Upper left, the first 2048 points, 
upper right the first 4096, down left, the first 6144 and down 
right, the first 8192. The units in all axes are given as multiples 
of 7T 



The frequencies covered by the power spectrum are 

Si = fc il j = 0-M 



(6) 



In the present work we used the so called "Hanning" window. 
More details on the calculation of the power spectrum can be 
found in the book by Press et. al. (1992). 

As a rule "mixed" methods are expected to perform bet- 
ter in distinguishing between ordered and chaotic trajectories. 
The reason for that is that time series that are constructed from 
the geodesic divergence of nearby orbits contain all the various 
characteristic frequencies that locally affect the motion in the 
"proper" ratio, i.e. the frequencies corresponding to the differ- 
ent directions (degrees of freedom) are properly weighted. In 
particular now, as far as our method is concerned, the Power 
Spectrum Of Divergences (PSOD) of an ordered trajectory is 
expected to posses only a few "spikes" at specific frequencies. 
The number of harmonics however depends on the system un- 
der consideration as well as on the values of its "controll" pa- 
rameters (see below). In contrast, the PSOD of a chaotic tra- 
jectory should appear continuous, due to the random nature of 
the qt time series. However, the above considerations lie be- 
hind all methods based on time series analysis. What is really 
important, for the assessment of the new method with respect 
to the other ones appearing in the literature, is to evaluate (a) 
its independence from the number of degrees of freedom of 
the dynamical system, (b) its effectiveness with respect to the 
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Fig. 2. The PSOD for the three trajectories of the 2-D mapping 
and for N=1024 iterations. Upper left : mapl trajectory, upper 
right : map2 trajectory, down left : map3 trajectory and down 
right : map3 trajectory but now with N=256 



object of test (single trajectories or distributions of initial con- 
ditions covering wide phase space regions) (c) its sensitivity, 
i.e. the minimum length of the time series, necessary to distin- 
guish a sticky chaotic trajectory from an ordered one and (d) 
its ability to produce a well-defined measure of chaos (as is the 
LCN). 

3. Evaluation of the method 

We proceed in the assessment of the method using three dif- 
ferent dynamical systems, namely a 2-D mapping as well as 
a 2-D and a 3-D Hamiltonian system. In each one we evalu- 
ated the nature of a considerable number of trajectories. In the 
following subsections we present only three or four trajecto- 
ries per system, which we think that are typical examples of 
the three different classes of trajectories, i.e. ordered, clearly 
chaotic and sticky. The amplitudes of the PSOD, in all figures, 
are normalized so that the highest has the value one, while the 
frequency is given in cycles per time unit. 

3.1. 2-D mapping 

We first test the method in the simple 2-D mapping 

Jj_l_i = Ji + k cos(26*i) mod(27r) 

i+ i = 0i + J i+X mod(27r) (7) 

where the stochasticity parameter k is taken equal to 0.7. 

We present here the results of three trajectories, an ordered 
one (mapl) starting at Jq = tt, 9q = 1.4 tt, a stochastic one 
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Fig. 3. The evolution of x(t) for the sticky trajectory in the 2-D 
mapping 

(map2) starting at Jo = 7r, 9q = 1.5 tt and a sticky one (map3) 
starting at Jo = tt, 9o = 1.538 tt. The initial conditions of 
the third one place it very close to the boundary between the 
ordered region, surrounding the stable point Jo = tt, 9o — 
57r/4, and the chaotic sea. Figure [l] shows the consequents (9, 
J) of the "sticky" trajectory at various times (iterations). As we 
can see, for up to 2 048 iterations the trajectory behaves like an 
ordered one. Around i = 4 000 it starts to present some signs 
of irregularity and finally, after i — 6 000, the chaotic nature 
of the trajectory becomes evident. Using our method on these 
three trajectories, with the length of the time series being 
N=1024 iterations and defining At = 1, we obtain the spectra 
shown in Fig. |[ 

The upper left frame of Fig. ^ is the spectrum of the or- 
dered trajectory. The spectrum consists of some basic frequen- 
cies while the "noise" is at a very low level. On the contrary, 
the spectrum of the chaotic trajectory, in the upper right frame, 
covers the whole frequency range with comparable amplitudes, 
i.e. a clearly continuous spectrum. Looking at the spectrum of 
the "sticky" trajectory (lower left frame), we see a pattern al- 
most the same as that of the chaotic one. Some high-amplitude 
spikes are evident but the continuum noise level is again very 
high for the whole frequency range, denoting the stochastic na- 
ture of the trajectory. Even with much less iterations (N=256 - 
lower right frame of Fig. |J) we get the same result. The spec- 
trum is less dense, since it has a smaller number of frequencies 
than before (see eq.(Q)), but the basic features are the same. 
Note that for the LCN, which is the limit of the function 

1 N 
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Fig. 4. The PSOD for the three trajectories in the 2-D mapping 
shown in Fig. |[ the peaks being sorted in descending order 
of amplitude. The first three frames are for N=1024 iterations, 
while the fouth ( down right ) is map3 trajectory with N=256 
iterations. 



as N — > oo (see Fig. |3|), or the plot of the consequents of the 
mapping (Fig. [I]), we need much more than 1024 iterations in 
order to decide whether the trajectory is chaotic or not. 

3.2. An improved criterion 

The above presented results show that it is, indeed, worth to 
consider the new method as a useful tool in the assessment 
of the nature, ordered or chaotic, of a trajectory. However the 
method, as it is, does not entail a clear and easy to apply crite- 
rion for the classification of a trajectory as ordered or chaotic. 
Here we try to improve somehow the presentation of the results 
of our method, in order to propose such a criterion. Note that 
this new criterion is similar, graphically, to the criterion of the 
FLI proposed by Froeschle et al (1997). 

If the peaks appearing in the PSOD are plotted in descend- 
ing order of amplitude, we have a graphical representation of 
how many strong frequencies the spectrum possesses. Figure Q 
shows this representation of the PSOD for the three trajectories 
studied in the previous sub-section. Ordered trajectories have 
only a few high amplitude frequencies and the background is 
formed by peaks whose amplitudes are more than four orders 
of magnitude smaller than that of the basic frequency. On the 
other hand the stochastic trajectories possess only a few high- 
amplitude frequencies and the largest part of the spectrum con- 
sists of a "continuum" of frequencies with also considerable 
amplitude. In this respect one may chose to represent the results 
by a signle number (e.g. the number of peaks up to a certain am- 
plitude), provided that a certain "threshold" for the noise level 
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Fig. 5. The PSOD for the regular (left) and the chaotic (right) 
trajectories in the 2-D mapping with N=256 (top), N=4096 
(middle) and N=65 536 (bottom). 



is chosen (10 in Fig. 4). This value will, however, depend on 
the system at study. 

3.3. Noise level 

Any method of analysis of finite-sample time series is bound 
to suffer from noise. Thus, even in the PSOD of a regular orbit 
a certain noise level is expected. This is mainly due to "leak- 
age" of power from the frequency lobes, a side effect of the 
calculation of the power spectrum using Discrete Fast Fourier 
Transform methods. When calculating the power spectrum of 
a "monochromatic" signal, the power contained in its basic 
frequency "leaks" into neighboring frequencies. The leakage 
depends on the windowing function used. When a signal pos- 
sesses two closely spaced frequencies, all frequencies in be- 
tween these two will also gain considerable amplitudes, due to 
this phenomenon. If the PSOD of a regular orbit has a large 
number of basic (strong) frequencies, then this effect can lead 
us to falsely identify it as chaotic. The problem can be tackled 
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Fig. 6. The sorted PSOD for the regular (left) and the chaotic 
(right) trajectories in the 2-D mapping for different iteration 
number, N. Lines from top to bottom correspond to N=256, 
1024, 4096, 16 384 and 65 536. 

at the cost of taking more points in the sample. In this way, the 
number of frequencies appearing in the spectrum is increased 
but the frequency lobes become thinner. Thus, lobe overlapping 
is reduced and the noise level drops. For chaotic orbits on the 
other hand, the observed noisy pattern is an inherent property 
of the spectrum and by increasing the number of points one 
cannot alter the picture. 

The above can be seen in Fig. |[ which shows the PSOD 
of a regular (left column) and a chaotic trajectory (right col- 
umn) for three different values of N, i.e. N=256 (top), N=4096 
(middle) and N=65 536 (bottom). As we see, the noise level in 
the case of the regular orbit drops significantly when N is in- 
creased, while in the case of the chaotic orbit it remains more 
or less unchanged. In the N=65 536 case the noise level for the 
ordered trajectory drops below 10 -20 , becoming comparable 
to the accuracy of the FFT calculation (double precision). 

This phenomenon is better seen if we use the amplitude- 
sorted PSOD. Fig. H shows the amplitude-sorted PSOD of the 
regular (left) and chaotic (right) trajectory with N=256, 1024, 
4096, 16 384 and 65 536. While the noise level in the PSOD 
of the ordered trajectory is reduced when N is increased, it re- 
mains the same for the chaotic one. 

3.4. 2-D Hamiltonian system 

We now proceed to test the method in a 2-D Hamiltonian 
system. We selected the Hamiltonian used by Caranicolas & 
Vozikis (1987), 

H=^(p 2 x +p 2 y )+x 4 + y 4 + 2ax 2 y 2 = h (9) 

where the parameter a is taken equal to 2.0 and the energy con- 
stant h = 1.0. 

Again we study three trajectories, one ordered {ord-2) start- 
ing at .t(0) = 0.5, one clearly stochastic (ch-2) starting at 
x(0) = 0.001 and one "sticky" (st-2) which starts at x(Q) = 
0.1326. All three trajectories have also y(0) = and p x (Q) = 
0, while their p lJ (0) is given by the energy integral. 

Using the PSOD with a renormalization time-step equal to 
At = 1.5, we find that the spectra of the three test trajecto- 
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Fig. 7. The evolution of x(t) of the four test trajectories in the 
2-D Hamiltonian model in the case a = 2. The initial positions 
of the trajectories are given in the text. Note that the \(t) curves 
of the two ordered trajectories have been lowered by one (ord- 
2) and half (ord-2a) orders of magnitude in order to be clearly 
distinguished from the sticky trajectory 




Fig. 8. The PSOD, with N=2048, of the four test trajectories in 
the 2-D Hamiltonian model for a = 2. Upper left: the ordered 
ord-2, upper right: the chaotic ch-2, lower left: the sticky st-2 
and lower right: the ordered ord-2a. 
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Fig. 9. The surface of section x, x of st-4 in the 2-D Hamilto- 
nian model. Upper left t — 3 000, upper right t = 8 000, down 
left t = 15 000, down right t = 20 000 



ries present the same properties as those of the corresponding 
cases of the mapping. In Fig. |8] we see the PSOD of the three 
test trajectories. The difference between the spectrum of the 
chaotic trajectory ch-2 and that of the ordered trajectory ord-2 
is again obvious. Moreover, the "sticky" trajectory st-2 has a 
spectrum similar to that of ch-2. Note that we have used only 
2048 points (corresponding to t — 3072) while, if we look at 
the evolution of the x(t) (Fig. ^), the trajectory looks ordered 
for a time up to t ps 7 000. The fourth frame (lower right) in 
Fig. H corresponds to another regular orbit (ord-2a) that starts 
at x — 0.133, i.e. very close to the sticky trajectory. We see 
that, although the two trajectories start very close to each other 
and, at least for the first 3 000 times steps, span approximately 
the same phase-space region, their PSODs are completely dif- 
ferent, clearly revealing the nature of each case. 

We decided to test also the case where the parameter a 
is taken equal to 4.0. As Caranicolas & Vozikis (1987) have 
shown, the surface of section of this case has a completely dif- 
ferent topology from that of the a = 2 case. The equipotential 
curves have negative curvature along the y — ±x lines. This 
affects mainly the loop orbits which appear "squared". 

Again we study four trajectories, one ordered (ord-4) start- 
ing at x(Q) = 0.3, one clearly chaotic (ch-4) starting at x(0) = 
0.7, one "sticky" (st-4) which starts at x(Q) = 0.36282 and 
one ordered (ord-4a) starting at x(0) — 0.362 very near to 
the sticky one. The surface of section plot for the "sticky" tra- 
jectory is shown in Fig. ^ at various times. Figures [K]and|n] 
present the PSOD and the amplitude-sorted PSOD for these 
four orbits (N=4096). An important characteristic seen in these 
two figures is the presence of a high number of medium- 



Fig. 10. The PSOD, with N=4096, of the four test trajectories in 
the 2-D Hamiltonian model for a = 4. Upper left: the ordered 
ord-4, upper right: the chaotic ch-4 and lower left: the sticky 
st-4 and lower right: the ordered ord-4a. 
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Fig. 11. The sorted PSOD of Fig. |To| 



amplitude frequencies in the spectra of the regular orbits. How- 
ever, a distinction between regular and chaotic orbits can still 
be made. Due to frequency overlapping, discussed in sec- 



tion 3.3, individual frequencies cannot be distinguished at an 
amplitude level smaller than 10~ 6 . This makes it very diffi- 
cult to identify a sticky orbit with an amplitude level around 
10~ 6 — 10~ 7 . Nevertheless the noise level of the regular or- 
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Fig. 12. The evolution of x{t) of the four test trajectories in the 
3-D Hamiltonian model. The initial positions of the trajectories 
are given in the text. 



bits will be supressed if we take more points, while for a sticky 
chaotic orbit it will remain more or less the same. 

3.5. 3-D Hamiltonian system 

We apply our method to the model Hamiltonian used by Mag- 
nenat (1982), Contopoulos & Barbanis (1989), Barbanis and 
Contopoulos (1995), Barbanis (1996), Varvoglis et al. (1997), 
Barbanis et. al. (1999) and Tsiganis et al. (2000a) 

H = \ (Pi + ?l + Pi) 

+ ^(Ax 2 + By 2 + Cz 2 ) - exz 2 -rjyz 2 = h (10) 

where the parameters are taken as A = 0.9, B = 0.4, C = 
0.225, e = 0.560 and r\ = 0.20 and the energy level is h = 
0.00765 

We again test three trajectories, an ordered one (regOl) 
starting at x — 0.01, y — 0.032, a chaotic one (ch-01) starting 
at x = 0, y — and a sticky one (st-01) with initial conditions 
x = 0.01725, y = 0.032, while z, p x , p y were taken equal 
to and p z is calculated from the energy integral. We use the 
variables x, y, z instead of x, y, z in order to be consistent with 
the previous publications. The barred variables are defined as 
x = \f~Ax, y — \f~By and z = \[Cz. 

In 3-D one cannot visualize a surface of section plot, in 
order to check whether a particular trajectory is ordered or 
chaotic. Therefore, if one is using the traditional tools, he has 
to rely on the calculation of LCNs. It should be pointed out 
that a positive LCN is a proof that the trajectory under study is 
chaotic, while a monotonically decreasing value of is not 
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Fig. 13. The PSOD of the four test trajectories in the 3-D 
Hamiltonian model using 512 points, i.e. 7680 time steps. Up- 
per left the ordered regOl, upper right the chaotic ch-01, lower 
left the sticky st-01 and lower right the ordered regOO. 

a proof of order, since this behavior could very well originate 
from "stickiness". That is why we decided to test one more tra- 
jectory (regOO) for which we can be almost certain that it is 
ordered, as it has the same initial position as regOl but it be- 
longs to an almost integrable case of the model Hamiltonian, 
e = 0.01 and?? = 0.01. 

Figure |l2| shows the calculation of the x(t) for the four 
trajectories. We can clearly see that x(t) of st-01 is decreasing 
up to t = 7 10 4 and then begins to saturate to a non-zero LCN 
value. The stochasticity is even more evident after t = 2.2 10 5 , 
where we have a "jump" to a higher LCN value. 

We calculate the PSOD using At = 15, a value approx- 
imately equal to the time interval between two consecutive 
crossings of the x — y plane by the trajectory. Figure [l3| shows 
the PSOD of the four trajectories using 512 points, i.e. for 
t = 7 680. Note once again that we can decide that the "sticky" 
trajectory st-01 is actually stochastic well in advance of the 
LCN method. The LCN shows the stochastic behaviour of the 
trajectory only after t — 70 000, while with the PSOD we need 
a modest t = 7 680. 

4. Comparison with other methods 

As we already mentioned in the introduction, three of the most 
recent methods for distinguishing ordered from chaotic trajec- 
tories are the Fast Lyapunov Indicators (FLI) (Froeschle et al. 
1997), the "spectra" of stretching numbers and/or twist angles 
(Froeschle et al. 1993, Voglis & Contopoulos 1994, Contopou- 
los & Voglis 1997) and the "spectral distance" (DSD) (Voglis 
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Fig. 14. The PSOD of the four test trajectories in the 3-D 
Hamiltonian model shown in Fig. [13], the peaks being sorted 
in descending order of amplitude. 




Fig. 15. The PSOD of the two trajectories used by Froeschle et 
al. (1997) to test the FLI method. In the left frame is the PSOD 
of the chaotic trajectory (x(0) = 0.001, y(0) = 0.001) and in 
the right frame is the PSOD of the ordered trajectory (x(Q) = 1, 

y(0) = 0). 



et al. 1998, 1999). Out of these three methods the fastest ones 
are the FLI and the DSD methods. In this section we compare 
the PSOD with the FLI and the DSD methods by applying it to 
the test models presented in the above mentioned papers. A dis- 
cussion concerning the spectra of stretching numbers follows. 

4.1. Single trajectories 
4.1.1. FLI 

In the paper by Froeschle et al. (1997) the authors tested the 
FLI method on two trajectories of the standard map 



Fig. 16. The PSOD (upper frames) and the sorted PSOD (lower 
frames) of the two test trajectories (left frames : regular A2, 
right frames : chaotic A3) in the 4-D mapping for N=4096. 

with k = 0.3; one stochastic starting at x = 0.001, y = 0.001 
and one ordered starting at x = 1, y = 0. They found that for 
the stochastic trajectory the FLI's drop very quickly down to 
10~ 20 (Fig. 2 in their paper - 200 iterations). On the contrary, 
the function x(t) levels only after about 10000 iterations. For 
the ordered trajectory the FLI's are slowly decreasing, follow- 
ing x(t)- I n Fig- H] we present the PSOD's of the stochastic 
(left frame) and the ordered (right frame) trajectory, calculated 
using only 256 iterations. It is obvious that the two spectra 
clearly differentiate between the two types of trajectory. 

4.1.2. Spectral distance (the 4D map) 

In a recent paper Voglis et al. (1999) proposed, as a tool for the 
distinction between chaotic and regular orbits in 4D maps, the 
use of the "spectral distance" D 2 . The method is based on the 
property that the "spectrum" of stretching numbers (as well as 
that of the helicity angles) of a chaotic trajectory is indepen- 
dent of the initial orientation of the deviation vector, while the 
spectrum of a regular trajectory is not. 

The "spectral distance" D 2 is a norm defined as 



D 2 =J2lSl(q)-S 2 (q)} 



(12) 



where the summation is for all q's and Si, S2 are two spectra of 
the same orbit but with two different initial deviation vectors. 

Voglis et al. (1999) applied their method to a 4-D mapping 
consisting of two coupled 2-D standard maps, i.e. 



Xi+i = Xi + ksm{xi + y{) mod(27r) 
Vi+l = Xi + y { mod(27r) 



(11) 
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(13) 
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where the Xi's are defined in the interval [0,1) (i.e.modl). 

We chose to test our method upon the two most interest- 
ing cases shown in Voglis et al. (1999), namely trajectories A2 
and A3, in their notation. The A2 case has initial conditions 
(xi, X2 1 X3, X4)=(0. 55,0. 1,0.62,0.2), (3 = 0.1 and is a regular 
orbit, while the A3 case has the same initial x t 's, (3 — 0.3051 
and is a chaotic orbit but with a very low value of LCN (around 
410" 7 ). 

Fig.|l6| shows the PSOD of the two test trajectories. The 
left panel corresponds to the regular orbit (A2) and the right 
panel corresponds to the chaotic one (A3). The spectra were 
calculated using 4096 iterations for each orbit. The distinction 
between the regular and the chaotic orbits is apparent. Note 
that our method gave the same result as the D 2 method with 
almost the same computational effort. Of course, both methods 
are much faster than the traditional LCN method. 

4.1.3. "Spectrum" of stretching numbers 

The "spectrum" of stretching numbers, S(q), (Froeschle et 
al. 1993, Voglis & Contopoulos 1994, Contopoulos & Voglis 
1997, Dvorak et al. 1998) is a method also based on the diver- 
gence of nearby trajectories. It consists of the calculation of the 
probability density of (eq.(|TJ)), i.e. 



S(q) = 



AN(q) 
Ndq 



(14) 



where AN(q) is the number of q/, with values between q and 
q + dq. A quasi-periodic trajectory, which lies very close to a 
periodic trajectory, has a "U" shaped distribution of stretching 
numbers which is also symmetric around q = (Contopoulos 
et al. 1997). As we move away from the periodic trajectory, 
this symmetry is destroyed (Caranicolas & Vozikis 1999) and 
the spectrum starts to develop a greater number of maxima. On 
the other hand if the trajectory is chaotic, the spectra have dif- 
ferent shapes and are not symmetric at all. It should be pointed 
out, however, that, in order to obtain a well defined spectrum, 
one needs to account for a large number of iterations (typically 
N = 10 5 or more). Therefore, we did not attempt to compare 
this method to our own. 

4.2. Sets of trajectories 

In order to circumvent the problem of calculating a trajectory 
for long times, Contopoulos and Voglis (1997) proposed the 
use of the average value of q, 



<q>- 



1 



JV 

\ " 



N ^ 



Ik 



(15) 



If we keep N small, we can scan a wide area of initial condi- 
tions and map its dynamical behavior. Trajectories in chaotic 
domains will have < q > scattered around the value of the 




Fig. 17. The logarithm of the average amplitude of PSOD in a 
cross section of 1000 trajectories near the 1/6 resonance in the 
standard map of eq.(6) 



LCN of this domain, while trajectories in the ordered domain 
will have < q > near zero. The method is very fast in distin- 
guishing between ordered and chaotic domains. Fig. 9 of Con- 
topoulos & Voglis (1997) is a very good example of the results 
obtained by this method with only 10 iterations. 

However, although the method is very good in scanning 
wide areas of phase space for locating islands of order, it cannot 
give reliable results with so few iterations for a particular tra- 
jectory. In the case of stochastic trajectories, < q > for N = 10 
varies so much, that it may yield a number as small as the one 
given for ordered trajectories. The situation is even worse in 
the case of sticky trajectories (i.e. in the borders of islands). In 
order to decide on the character of such trajectories one needs 
considerably longer calculations. 

Froeschle & Lega (1998) tested the FLI method along with 
the method of twist angles (Contopoulos & Voglis 1997), the 
frequency map analysis method (Laskar et al. 1992, Laskar 
1993) and the sup-map method (Laskar 1990, Froeschle & 
Lega 1996). Figs. 8a-d of their paper shows the results of the 
four methods on a cross section of 1 000 trajectories near the 
hyperbolic point of the 1/6 resonance of the standard map (Eq. 
11) with k = —1.3. For the FLI method they used 2000 iter- 
ations while for the other three methods 20000 iterations. In 
order to produce unambiguous results, the trajectories in this 
test are classified as ordered or chaotic by an appropriately se- 
lected number/indicator. For the FLI method the authors have 
used as an indicator the time necessary for the FLI to reach a 
value lower than 10~ 10 . 

In order to compare our method with these results we need 
also an one-number indicator derived from the PSOD. As such 
we have selected here the average value of < P/P max >■ The 
averaging is performed not over all values but by ignoring the 
highest l/6th and the lowest l/6th of the amplitudes, the former 
being probably due to periodicities in sticky trajectories and the 
latter probably coming from numerical errors in the integration 
and the calculation of the power spectrum. Fig. [l7| shows the 
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Fig. 18. The evolution of x(f) of the trajectories of four aster- 
oids. Upper left, an ordered trajectory with initial a = 3, e = 0; 
Upper right a highly stochastic trajectory with a = 4, e = 0.2; 
Lower left, a stochastic trajectory with stickiness and with ini- 
tial a = 3.64, e = 0.08; Lower right, another sticky orbit with 
initial a = 3.63, e = 0.08 



same cross section with Fig. 8a-d of Froeschle & Lega (1998), 
using the above indicator taken after 8 192 iterations. As we can 
see it gives essentially the same information as the other four 
methods. Note that the y-axis in Fig. [l7| is inverted for easier 
comparison with Fig. 8a-d of Froeschle & Lega (1998). 

5. Application to asteroidal motion 

As an application of our method to a problem of physical im- 
portance, we shall use it in order to assess the chaotic or not na- 
ture of asteroidal trajectories. We use a simplified model of the 
solar system, namely the planar restricted three body problem 
where the Sun and Jupiter move on elliptic trajectories around 
their center of mass and an asteroid of infinitesimal mass moves 
in the gravitational field of the two bodies. We calculate the 
PSOD and the LCN using as renormalization time, At, the pe- 
riod of Jupiter, i.e. At = Tj ~ 11.86 years. 

Fig. |l8| shows the x{t) of the f° ur trajectories tested. The 
upper left frame belongs to an ordered trajectory starting with 
semi-major axis a = 3 and eccentricity e = 0, the upper right 
to a stochastic trajectory with initial elements a = 4, e = 0.2. 
The other two frames correspond to trajectories with initial a = 
3.64, e = 0.08 (lower left) and a = 3.63, e = 0.08 (lower 
right) representing "clones" of the asteroid 522 - Helga, which 
is a well-known example of "stable chaos" (Milani & Nobili 
1992). 




Fig. 19. The PSOD of the four asteroidal trajectories with 
N=512, i.e. 6062 years 




Fig. 20. The PSOD, with N=8192, of the ordered asteroidal 
trajectory (left) and the chaotic Helga clone with a — 3.63, 
e = 0.08 (right) 



The PSOD of the four trajectories (N=512) is shown in 
Fig. [l9|. As we can see, a few Jovian periods are sufficient to 
decide if the motion of a particular asteroid is stochastic or or- 
dered. In the case of the Helga clone with a = 3.64 (lower left 
of Fig. [1^) the PSOD shows clearly a chaotic nature after only 
512 Jovian periods i.e. 6 072 years, while even a rough calcula- 
tion of the LCN needs at least 10 5 years (see lower right frame 
of Fig. [jj|). The PSOD of the a = 3.63 Helga clone is rather 
peculiar. Although it differs from that of the ordered orbit (up- 
per left), it is not clearly chaotic, unlike the PSOD of the other 
Helga clone (lower left). Nevertheless, if we take more points 
in our time series the chaotic nature of the orbit becomes ap- 
parent, as we can see in Fig. However, the spectrum can 
still be described as having a strong quasi-periodic component, 
something which is related to the peculiar dynamical nature of 
this orbit as discussed in Tsiganis et al. (2000b). 
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6. Conclusions - Discussion 

In the present paper we propose an alternative tool, which we 
call PSOD, for the characterization of the chaotic or not nature 
of trajectories in conservative dynamical systems. The method 
is based on the frequency analysis of a time series, constructed 
by successive records of the amplitude of the deviation vec- 
tor of nearby trajectories. As discussed in Section 2, such a 
"mixed"[] method is expected to have certain advantages. The 
reason is that the power spectrum of such a time series will con- 
tain all the characteristic frequencies of the motion in a prop- 
erly "weighted" ratio. 

The basic characteristic of the PSOD, seen in all three test 
cases (a 2-D mapping, a 2-D and a 3-D Hamiltonian system), 
is that 

- for ordered trajectories the spetrum possesses only a few 
high-amplitude peaks, the exact number of which depends 
not only on the system but also on the particular orbit. Of 
course a small amplitude noise level, due to the numerical 
procedure, is superimposed on the spectrum, which dimin- 
ishes as the length of the time series is increased; 

- for chaotic trajectories the spectrum has a noisy pattern. 
For weakly chaotic orbits a few high-amplitude peaks are 
also present. Increasing the length of the time series, the 
spectrum tends to a white noise spectrum which remains 
practically unchanged for any (large enough) number of 
points (see Figs. 5 and 6). 

Like most methods existing in the literature, it seems that 
the method performs better for maps than for flows. However, 
the results found for the three Hamiltonian flows tested (includ- 
ing the three-boody problem) show that the method can be ap- 
plied to any system, no matter how many degrees of freedom 
are involved. We believe that the results may be significantly 
improved for flows, provided that proper analysis concerning 
the renormalization time is conducted. The difference between 
maps and flows is that, in the former case, isochronous records 
of dynamical quantities also mean studying the system on a sur- 
face of section. This is not the case for flows and, thus, a more 
refined analysis on how to select a proper renormalization time 
has to be made. 

We have shown that the sensitivity of the PSOD in testing 
single trajectories in maps is comparable to the FLI and the 
DSD methods. This also makes the PSOD an efficient tool for 
scanning wide areas of the phase space (see Section 4). For 
such purposes one would like to have an one-number indica- 
tor for measuring chaos. The only uniquely defined measure of 
chaos is of course the LCN. Any other indicator should be in a 
one-to-one correspondence with the LCN, in order to give the 
same information. There is no guarantee that such an indicator 
can be based on the frequency content of the PSOD, as chaotic 
orbits with similar LCNs may have a different frequency dis- 
tribution. Further analysis of the properties of the PSOD has to 



be performed in order to decide whether such an indicator can 
be found. 
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